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Abstract. We derive a numerical method for Darcy flow, hence also 
for Poisson's equation in first order form, based on discrete exterior cal- 
culus (DEC). Exterior calculus is a generalization of vector calculus to 
smooth manifolds and DEC is one of its discretizations on simplicial 
complexes such as triangle and tetrahedral meshes. We start by rewrit- 
ing the governing equations of Darcy flow using the language of exterior 
calculus. This yields a formulation in terms of flux differential form and 
pressure. The numerical method is then derived by using the framework 
provided by DEC for discretizing difl'erential forms and operators that 
act on forms. We also develop a discretization for spatially dependent 
Hodge star that varies with the permeability of the medium. This also 
allows us to address discontinuous permeability. The matrix represen- 
tation for our discrete non-homogeneous Hodge star is diagonal, with 
positive diagonal entries. The resulting linear system of equations for 
flux and pressure are saddle type, with a diagonal matrix as the top left 
block. The performance of the proposed numerical method is illustrated 
on many standard test problems. These include patch tests in two and 
three dimensions, comparison with analytically known solution in two 
dimensions, layered medium with alternating permeability values, and a 
test with a change in permeability along the flow direction. A short in- 
troduction to the relevant parts of smooth and discrete exterior calculus 
is included in this paper. Wo also include a discussion of the boundary 
condition in terms of exterior calculus. 



1. Introduction 

We have discretized the equations of Darcy flow and obtained a numer- 
ical method on staggered mesh pairs with fluxes and pressures being the 
primary variables. The numerical method was obtained by using discrete 
exterior calculus (DEC) [22, 23, 28]. Exterior calculus generalizes vector cal- 
culus to higher dimensions and to smooth manifolds [1] and DEC is one of its 
discretizations. This discretization yields numerical methods for solving par- 
tial differential equations (PDEs) on simplicial complexes, such as triangle, 
tetrahedral or higher dimensional simplicial meshes. A recent implemen- 
tation of DEC is described in [7] . DEC is related to many discretizations 
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of exterior calculus that have been popular in or are currently being pur- 
sued in numerical analysis. Others include finite element exterior calculus 
[4], mimetic and compatible discretizations of PDEs [10, 31], and the covol- 
ume method [40]. See [3] for a collection of recent papers in these fields. 
Those parts of exterior calculus and DEC that are relevant to this paper are 
summarized in Section 2. 

The equations of Darcy flow model the flow of a viscous incompressible 
fluid in a porous medium. The equations consist of Darcy's law (which ex- 
presses force balance), the continuity equation, and the boundary condition. 
The first two form a very simple pair of equations, being Poisson's equation 
in first order form. In this case Darcy flow can be rewritten as Poisson's 
equation with pressure as the unknown. In applications however, it is often 
the velocity field or flux that is of primary interest and there is a loss of 
smoothness and accuracy in going from pressure to velocity. In flrst order 
form (i.e., when Darcy law and the continuity equations are not combined 
into a single equation) Darcy flow equations can be discretized using mixed 
flnite element or volume methods. 

In mixed methods, velocity (or flux) along with pressure are taken to be 
the primary variables and this can yield more accurate results compared 
to the pressure-only formulation. Mixed formulations require careful choice 
of spaces for velocity and pressure since not all combinations yield stable 
methods. For example, the use of continuous piecewise linear representation 
for both velocity and pressure results in an unstable method [16, 18]. An 
example of such unstable behavior is shown in Figure 1. Many fixes for such 
instability are presented in the literature. For example, one can use Raviart- 
Thomas (RT) elements [42], Nedelec elements [39], Brezzi-Douglas-Marini 
(BDM) elements [17], or a variety of other finite element spaces summa- 
rized in [19]. Many of these spaces that have yielded stable discretizations 
of Darcy fiow and other problems have been unified under the umbrella of 
finite element exterior calculus [4] . Stabilized mixed finite elements methods 
have also been developed for Darcy flow [9, 30, 35, 37, 38]. Finite volume 
[2] and covolumc methods [20, 21, 40] arc yet another approach for solving 
the equations of Darcy fiow. A monotone, locally conservative finite vol- 
ume scheme for diflFusion equation is described in [34] and a mimetic finite 
difference scheme is in [33]. 

The primary variables in our method are area or volume fiux (depending 
on whether the problem is 2D or 3D), and pressure. The fluxes are placed 
on edges in triangle meshes or triangles in tetrahedral meshes. This leads to 
pressures being placed at circumcenters of top dimensional simplices. Thus 
the pressure can be considered constant in each simplex. Our method is 
related to the methods for diffusion described in [41] and to the covolume 
method applied to Darcy flow in [20, 21] but our emphasis is on the rela- 
tionships between the numerical method and exterior calculus. In addition 
we treat discontinuous permeability explicitly. The treatment of flux and 
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Figure 1. Mixed finite element method can be used to solve Darcy 
flow and it is well-known that care is needed in selecting the under- 
lying finite element spaces. For example, equal order interpolation 
for both velocity and pressure is unstable. Here we show results 
for such a choice using piecewise linear finite elements. The correct 
solution is the constant vector field (1,0). The fix for such unsta- 
ble behavior is well-known in the finite elements literature [16, 18]. 
In this paper we provide a related but different numerical method 
based on discrete exterior calculus. 



pressure in our method is similar to that for Navier-Stokes equations in [26], 
which is also based on DEC. 



Our results: Our discretization space is similar to the lowest order Raviart- 
Thomas elements. However, the linear system matrix resulting from our dis- 
cretization is a saddle type matrix [8] with a diagonal matrix as the top left 
block. Our method enforces mass balance locally and globally and passes 
several standard numerical tests. In 2D we show that the method passes 
patch test involving constant velocity and linear pressure. We also compare 
our numerical solution with an analytical solution for a more general source 
term. We develop a diagonal discretization of spatially dependent Hodge 
star that depends on the permeability. This is used for the case of layered 
medium with alternating permeabilities, and for a discontinuous medium 
where we use different pairs of permeability under constant velocity con- 
ditions. We also show that our method passes patch tests in the 3D case. 
These numerical results are in Section 6 and the advantages of a DEC based 
approach are discussed in Section 7. A discussion of boundary conditions in 
terms of exterior calculus is in Sections 3 and 6.1. 
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2. Review of Discrete Exterior Calculus (DEC) 

In this section we briefly outline the relevant parts of smooth and discrete 
exterior calculus. We discuss only the operators that are relevant for Darcy 
flow. For more details on DEC see [22, 23, 28] and for details on exterior 
calculus see [1, 5]. In Sections 4 and 5 we describe our method for solving 
equations of Darcy flow so that it can be implemented without knowledge 
of DEC or exterior calculus. Thus, a reader unfamiliar with some of the 
terms used in this section should still be able to follow and implement the 
method. One useful characteristic of exterior calculus is that all objects and 
operators can be expressed in coordinate independent fashion. This aspect 
however is harder to explain in a few paragraphs. Instead, we give some 
examples using coordinates to describe the operators and objects of exterior 
calculus. 

2.1. Smooth exterior calculus. As mentioned earlier, exterior calculus 
generalizes vector calculus to smooth manifolds [1, 5] and it consists of op- 
erators on smooth general tensor fields defined on manifolds. A tensor field 
evaluated at a point is a multilinear function on the tangent space, mapping 
vector and covector arguments to M (the set of real numbers). Other ranges 
besides M are possible but not relevant here. Vector fields, symmetric tensor 
fields such as metrics and antisymmetric tensor fields are all examples of 
tensors. Antisymmetric tensors have been singled out in exterior calculus 
and are called diflFerential forms. A differential fc-form when evaluated at a 
point is an antisymmetric multilinear map on the tangent space that takes 
k vector arguments and produces a real number. It is an object that can 
be integrated on a A;-dimcnsional space. In exterior calculus, it only makes 
sense to integrate differential A:-forms on a fc-manifold. 

Let M be an n-dimensional orientable Riemannian manifold (a manifold 
with an inner product on the tangent space at each point), TM the tangent 
bundle (disjoint union of the tangent spaces at all points of M), X(M) the 
space of smooth vector fields and Q^{M) the space of differential A;-forms 
on M. 

Then exterior derivative is a map dfe : n''{M) n''+^{M) (sometimes 
written without the subscript) that raises the degree of a form, and the 
wedge product is a map or binary operator A : ri^(M) x r2'(M) — > r2'^+'(M) 
that combines differential forms. The most important property of d is that 
dfc+i o dfc = 0. These two operators are enough to describe a basis for 
differential forms on M. Taking M = (the standard three dimensional 
Euclidean space), with standard metric and coordinates x, y and z, a basis 
for r2^(R^), the space of 1-forms is {dx,dy,dz) and a basis for Jl^(R^) is 
{dx A dy, dx A dz, dy A dz). Let / be a scalar valued function on (i.e., a 
0-form) . Then its exterior derivative d / equals its differential df and is 
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For a 1-form a = aidx + a2dy + asdz, where ai are scalar valued functions, 
its exterior derivative is 

f da2 dai\ f das 9ai\ f das ^"2 \ , ^ , 

= ( - ^^"^^^^+1 5^ - 97 J 97 - -di ) ^^^^^ ■ 

The Hodge star is an isomorphism, * : il^(M) — r2"~'^(M). For with 
standard metric, the Hodge star satisfies the following properties: 

* 1 = dx A dy A dz ; 

*dx = dy A dz, *dy = —dx A dz, * dz = dx A dy ; 
*{dy A dz) = dx, *{dx A dz) = —dy, *{dx A dy) = dz ; 
*{dx Ady A dz) = 1 . 

For the equivalent properties are 

*l = dxAdy; 
*dx = dy, *dy = —dx ; 
*{dx A dy) = 1 . 
An important property of Hodge star is that for a /c-form a, 

(2.1) **a = 

Another operator relevant for Darcy flow is flat, which is a map b : 
X{M) ri^(M) that identifies vector fields and 1-forms via the metric. 
Consider with the standard inner product, and standard orthonormal 
basis. Then for a vector field V with components Vi, V2 and V3, we have 
V'" = Vidx + V2dy + Vsdz . If the inner product is not the standard one or the 
basis is not orthonormal then the relationship between a vector field and its 
flat in coordinates is more complicated. In M^, for a scalar function / and 
vector field V, some important relationships involving flat operator are: 

(2.2) (V/)^ = d/, {cmivf = *dV\ divF = *d*y^ 

Thus div, grad and curl can be defined in terms of exterior calculus opera- 
tors. Note that the operators d and A are metric independent and so they 
can be defined on a manifold without having to define a Riemannian met- 
ric. On the other hand, the operators * and b do require a metric for their 
definition. 

2.2. Primal and dual mesh. Discretizing exterior calculus involves decid- 
ing what should replace the smooth manifolds, differential forms and other 
tensor fields and operators that act on these. Recall that a simplicial com- 
plex K in M-^ is a collection of simplices in such that every face of a 
simplex of X is in X and such that the intersection of any 2 simplices of K 
is a face of each of them. The dimension n of the complex is the highest 
dimension of its simplices. Thus n < N where A'^ is the dimension of the 
embedding space. In DEC, the oriented Riemannian manifolds M of smooth 
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exterior calculus is replaced by the underlying space \K\ of an oriented man- 
ifold simplicial complex K as described in [28]. Briefly, these are simplicial 
complexes in which the neighborhood of every interior point is homcomor- 
phic to ("looks like") an open subset of and in which each simplex of 
dimension A;, for < fc < n — 1, is a face of some n-simplex. Thus a triangle 
with an edge sticking out from one vertex would not be admissible and nei- 
ther would a triangle mesh surface with a fin like triangle sticking out from 
an edge. In addition all the n-simplices must have the same orientation. See 
[28] for details. An oriented manifold simplicial complex will be also called 
a primal mesh. 

In addition to the primal mesh K a staggered cell complex ★ K associated 
with K and referred to as the dual mesh also plays a role in DEC. An 
example of a primal mesh, with some pieces of the dual mesh highlighted is 
shown in Figure 2. Usually the dual mesh is not explicitly stored. What are 
needed instead are lengths, areas or volumes of pieces of the dual mesh, the 
specific needed quantities depending on the application. The primal mesh is 
a simplicial complex, i.e., a triangle or tetrahcdral (or higher dimensional) 
mesh such as are used in finite volume or finite element methods. The only 
requirement on the primal mesh is that it be a Delaunay triangulation [25] . 
This is a requirement that most modern meshing programs automatically 
meet. Some examples are the Triangle [43] and tetgen [44] programs. 

In what follows, will be a primal fc-simplex, a fc-dimensional simplex in 
the primal mesh. The corresponding dual (n— A;)-cell, an (n— A;)-dimensional 
cell in the dual mesh will be denoted by *(7*^. We will use cc(o-'^) to mean 
the circumcenter of o"^, the unique point equidistant from all vertices of . 
The notation a < t will mean that simplex cj is a face of simplex r and 
T )^ a will mean that r contains a as one of its faces. 

To find the dual cell -ka^ of the primal simplex proceed as follows. 
Start from the circumcenter cc((T^). Traverse in straight lines, one by one, 
to the circumcenters cc((j'^+^) of all a^^^ >- and from those to the next 
dimension and so on all the way to the top dimension. Each path of these 
tr aver sals yields a simplex of dimension n — k. The union of all such simplices 
is the dual cell -ka^ of simplex a^. For example, in Figure 2, the dual of an 
internal edge is obtained by starting from its middle and traversing to the 
circumcenters of the two triangles containing it. The resulting two edges 
together form the dual of the edge. Note that it need not be a straight line 
if the two adjacent triangles do not lie in the same affine space. See [28] for 
more examples of primal-dual pairs. 

The primal and dual meshes of DEC are oriented. The primal mesh is 
oriented consistently at the top level. For example, either all the triangles 
in a triangle mesh must be oriented clockwise, or all of them must be ori- 
ented counter-clockwise. Similarly, all the tetrahedra in a tetrahedral mesh 
must be right-handed, or all must be left-handed. The lower dimensional 
simplices can be oriented arbitrarily, for example, using the dictionary order 
of the vertex numbers. The orientation of the dual cells is implied by the 
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Figure 2. A simplicial complex is subdivided using circumccn- 
tric subdivision and the dual cells are constructed from the sub- 
division. The new edges introduced by the subdivision are shown 
dotted. The dual cells shown are colored red. The dual cells are 
not explicitly stored, only their length, area, or volume is needed. 
The only requirement on the primal mesh is that it be a Delaunay 
triangulation, a requirement that is met by most modern mesh- 
ing software, such as Triangle [43] and tetgen [44]. Figure taken 



orientation of the corresponding primal simplices and of the top level primal 
simplices. For details see [28]. For triangle meshes, a vector along the primal 
edge followed by one along the dual edge should define the same orientation 
as that of the triangles. For tetrahedral meshes the orientation of a face 
followed by a vector along the dual edge should form the same handedness 
as that of the tetrahedron. 

2.3. Chains and boundary operator. Let be a finite simplicial com- 
plex. Recall from algebraic topology [36] that a k- chain on if is a function 
from the set of oriented A:-simplices of K to the set of integers Z. A A:-chain 
c has the property that that c{a) = — c(— cr), where —a is a with the op- 
posite orientation. Chains are added by adding their values. The space of 
/c-chains is a group and is denoted Ck{K). The group structure will not 
be important to us except for the fact that it will allow us to talk about 
homomorphisms - maps between groups that preserve the group structure. 
For a /c-simplex we will use or a to denote the simplex as well as the 
chain that takes the value 1 on the simplex and on all other /c-simplices 
in K. Such a chain is called an elementary k-chain. 

The boundary operator : Ck{K) — )■ Ck-i{K) is defined as a homomor- 
phism (that is, dk[a + h) = d^a + b) by defining it on oriented simplices 
[vq, . . . Vk]- It is defined by 



from [28]. 
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the hat indicating the missing vertex. For example if [fo,'Vi,i'2] is a tri- 
angle in M^, oriented counterclockwise, then the boundary of the corre- 
sponding elementary 2-chain is the sum of the 3 elementary 1-chains [ui, ^2]) 
-[vo,V2]{= [v2,vo]) and [vo,vi]. 

2.4. Cochains as discrete forms. As is usual in most discretizations of 
exterior calculus, in DEC discrete differential fe-forms are defined to be ele- 
ments of Hom(Cfc(-fC), M) . This is the space of real- valued homomorphisms 
on the space of fc-chains. This space is called the space of k- cochains or dis- 
crete differential k-forms. Thus given A;-chains a and b and a /c-cochain a, we 
have a{a-\-b) = a{a) +a{b) in which each term is a real number. The space 
of primal fc-cochains on a simplicial complex K will be denoted by C^{K; M) 
and the dual fc-cochains on the dual cell complex -kK hj C^{k K]M). We 
will shorten these notations to C^{K) and C''{~kK). On the underlying 
space \K\ of the simplicial complex K, we will be working with piecewise 
smooth differential forms. More precisely DEC starts with square integrable 
forms [24] on 1^1, denoted as L^Q''{\K\) and so that is the notation we will 
use for piecewise smooth forms. Discrete forms are created from (piece- 
wise) smooth forms with the de Rham map R : L'^Q,^{\K\) — )■ C^{K) or 
R : L^O'^diiTl) — )■ C^{'kK) depending on the context. See [24] for details 
on the de Rham map. For a form a, we will denote the evaluation of the 
cochain R{a) on a chain c as {R{a),c) and define it as J^a. Thus given a 
smooth fc-form a the de Rham map converts it into the fc-cochain J a with 
the slot for integration domain left empty. This cochain J a is ready to be 
applied to a A;-chain c to produce a number a. Note that we are implicitly 
assuming that the smooth quantities are defined on the simplicial mesh that 
is the discretization of the smooth manifold. For planar and spatial domains 
we consider in the examples in this paper, this condition is trivially true. 
For more general domains, like surfaces embedded in this restriction can 
be removed, but work on such generalizations, especially how it affects nu- 
merical solutions of PDEs is still in early stages. For some related ideas see 
[27, 48]. 

The cochain spaces C^{K) and C^{-kK) defined above are groups, but in 
addition they are also vector spaces. If there are A^^ simplices of dimension 
k m. K then the vector space dimension dim (C*^ (AT)) is A^^. Similarly, if 
there are ATj. cells of dimension k m. -kK then dim.{C^{-kK)^ is N^. Note 
that d\m.(C^{K)^ = dim(C"'~'^(* -fC)) since /c-simplices of K arc in one-to- 
one correspondence with (n— /c)-cells of ★AT. To define a basis for C^{K) as 
a vector space, the fc-simplices are first ordered in some way. For example 
in PyDEC [7] , the A;-simplices are ordered in dictionary order based on the 
vertex names of the simplex. In a triangle [fo, ^^i, ^^2]) the dictionary order 
for edges is [vcfi], [^'0,^2] and [^1,^2]- We would typically refer to these 
edges as ctq, a\ and a\, respectively. If ctq, . . . , ctn^. is an ordered list of the 
/c-simplices of K then we can define a basis {ctq, . . . , c^,, ) for C^{K) where 
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{a* , Uj) = 6ij, the Kronecker delta. That is, a* is the /c-cochain that is 1 on 
the elementary A;-chain of ai and on the other elementary fc-chains in K. 
In computations we represent elements of C^{K) and C^{-kK) as vectors of 
appropriate dimensions in these bases. 

2.5. Discrete exterior derivative and Hodge star. We now give def- 
initions of discrete exterior derivative and discrete Hodge star. These are 
given here without explanation as to why these are good choices for the dis- 
crete operators. Such an explanation can be found in [22, 28]. The discrete 
exterior derivative on the primal cochains will also be denoted as d (or dfc 
if the degree of the source space is to be specified) and is defined using the 
boundary operator in such a way that Stokes theorem is true by definition. 
For the exterior derivative on the dual cochains we will use the notation d* 
(or d^). For a fc-cochain and {p + l)-chain c^'^^, define 

(2.3) (dfca^c^+l):=(a^afe+lc'=+^). 

In the basis for C^{K) described in Section 2.4, the matrix representation 
of dfc, denoted Df^, is an N^j^i by A^^ matrix with entries 0, 1 or —1. Anal- 
ogous bases for C^{kK) using the dual cells yield matrix form of the dual 
discrete exterior derivative operator. The matrices of the primal and dual 
exterior derivative arc related. In fact the matrix form of d^, is ibl?^_^_^ 
where the sign depends on n and k. For Darcy flow the only relevant pairs 
of primal and dual exterior derivatives are dg and dn-i and the matrix form 
for dl is (-l)-^'^!- 

The discrete Hodge star can be thought of as an operator that "transfers 
information" between the primal and dual meshes. Given a A;-simplex a and 
a primal A;-cochain a, the discrete Hodge star of a (denoted * a) is a dual 
discrete (n — A;)-cochain defined by its value on the dual cell * cr by 

(2.4) -^(*a,*(j) := 7^(a, (j) . 

Here \a\ is the measure of a and |*(t| is the measure of the circumcentric dual 
cell corresponding to u, with measure of a 0-dimensional object being 1. See 
[22, 28] for details. We will sometimes use *fc to denote the Hodge star with 
C^{K) as its domain. Using the bases for C^{K) and C^~^{-kK) mentioned 
above, the matrix for *k is a diagonal Nk by matrix. It has been 
shown recently [29] that for a Delaunay mesh the diagonal entries are all 
positive. It is helpful to use to denote the inverse map although the 
inverse notation is usually not used for smooth Hodge star. To simulate the 
Hodge star property (2.1) on the discrete side we will always use 

(2.5) (-l)^("-^^fc\ 

whenever the inverse discrete Hodge star map is used. 

The various primal and dual cochain complexes in 2 dimensions are related 
via the discrete exterior derivative and discrete Hodge star operators in the 
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manner shown in the fohowing diagram: 



C^{K) C\K) C^{K) 



di 



(2.6) 



*o 



*i 



*2 



The corresponding diagram for 3 dimensions is: 

CQ{K) C\K) C^{K) C^K) 



(2.7) 



*o 



*i 



*2 



*3 



2.6. Interpolation of cochains. To go from cochains to piecewise smooth 
forms on a simpHcial complex K a well-known map is the Whitney map 
[11, 13, 49] denoted W. This map can be thought of as a way to interpolate 
numbers defined on edges, triangles and tetrahedra. Thus it goes in direction 
opposite of the one for the de Rham map which is a discretization map. For 
a complex K consisting of a triangle [i'0)^i)^2] and its faces, the Whitney 
map for 1-cochains is defined by extending by linearity, the following to all 
of C^K) 

(2.8) W{[vi,Vj]*) := m djij-Hj dm, 

for i ^ j and i,j G {0,1,2}. Here is the barycentric basis function 
corresponding to Vi, i.e. the affine function that is 1 on Vi and at other 
vertices. Recall that the 1-cochain [vi, Vj]* is 1 on the elementary 1-chain of 
the edge and at all other elementary 1-chains. For a tetralicdron 

[vo,vi,V2,V3] there are Whitney maps as above for interpolating the edge 
values and in addition there are Whitney maps for interpolating the triangle 
values, and these are defined by extension from 

(2.9) W{[vi,Vj,Vk]*) := 2{fii dfij Adfi^- jJ-j d^Ad^k+l^k d^ Adfij) . 

The piecewise smooth forms (smooth in each simplex) constructed using 
the Whitney map are also known as a Whitney forms. Whitney forms can 
be used to build a low order finite element exterior calculus [12, 13] and in 
computational electromagnetism the Whitney 1-form and 2-form arc also 
known as edge and face elements respectively [14]. Finite element exterior 
calculus has now been generalized to general polynomial differential forms 
[4]. We use the Whitney maps only for interpolating the differential forms 
so we can plot the corresponding vector field for visualization. The vector 
field corresponding to the Whitney 1-form in equation (2.8) is obtained by 
applying a sharp operator to get 
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k 


in M, 


V H — vp 


= -P9 








divf 


= 4> 


in M, 


V ■ h 




on dM, 



The vector field corresponding to the Whitey 2-form in equation (2.9) is 
2(//j V/Xj X V/Xfc - iij Vjii X Viik + Aife V//i X Viij) . 

3. Governing Equations 

We first present the governing equations of Darcy flow in the standard 

vector calculus notation, and then rewrite them using differential forms and 
vector fields (that is, in exterior calculus notation). This latter form is then 
discretized on a simplicial complex and its dual, which yields a numerical 
method for Darcy flow. 

Let M C be a bounded open domain, M its closure and dM := M\M 
its boundary, which is assumed to be piecewise smooth. In this paper n 
(which represents spatial dimensions) can be 2 or 3. Let v : M — > be the 
Darcy velocity [35] (units m^/(ms) = m/s for n = 2 or m'^/(m^s) = m/s 
for ra = 3) and let p : M ^ M be the pressure. The governing equations of 
Darcy flow can be written as 

(3.1) 

(3.2) 
(3.3) 

where A; > is the coefficient of permeability of the medium (units m^ for 
n = 3), /i > is the coefficient of (dynamic) viscosity of the fluid (units 
kg/(ms)), p > is the density of the fluid, g is the acceleration due to 
externally applied body force (i.e., pg is the body force density), ^ : M ^ M 
is the prescribed divergence of velocity, ^ : dM — )■ M is the prescribed normal 
component of the velocity across the boundary, and h is the unit outward 
normal vector to dM. For consistency cj) dM = Jqj^ dT where dV is the 
measure on dM. 

Equation (3.1) is Darcy's law, equation (3.2) is the continuity equation 
and equation (3.3) is the boundary condition. In the above equations, per- 
meability k is assumed to be a scalar constant. In Section 5 we will relax this 
constraint and allow /c to be a scalar valued function of space. The further 
generalization needed for modeling anisotropic permeability requires k to be 
a tensor, which is not addressed in this paper. To simplify the treatment, in 
the rest of this paper we will assume that there is no external force acting 
on the system, i.e., g = 0. 

The first step in the DEC formulation is to rewrite the governing equations 
(3.1) (3.3) in exterior calculus notation. As above, we first assume that 
the permeability is a scalar constant. We first apply the fiat operator 
to both sides of equation (3.1), use equation (2.2) for divergence, and then 
apply Hodge star to both sides of equation (3.2) to obtain (assuming g = 0) 

k 

(3.4) v' + - do p = in M , 

A* 
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(3.5) dn-ii*v^) = (l)U} inM, 

(3.6) *v^='>p'j ondM. 

Here u; = * 1 is a volume n-form on M and 7 is tlie volume (n — l)-form on 
the boundary dM and it is defined by requiring 

(3.7) liXi, . . . ,X„„i) = uj{h,Xi, . . . ,X„_i) , 

for all vector fields Xi, . . . , Xn-i on the boundary dM. Note that in going 
from equation (3.3) to (3.6) we have used the fact that {v ■ n)7 = . For 
an explanation of why this is true see [1, page 506]. The definition (3.7) of 
7 has implications on how the orientations affect the sign of the quantity 
-07 and this is explained using a concrete example in Section 6. The other 
quantities in the equations above are as in (3.1) - (3.3). Note that ^ and ^/j 
must satisfy 



M JdM 

by Stokes' theorem, which is analogous to the consistency condition stated 
earlier. 

Next we define a differential form which is the volumetric (or volume) flux 
in 3D (units m^/(ni^s) = m/s) or area flux in 2D (units m^/(ms) = m/s). 
This quantity will be denoted as /, which is an (n — l)-form, defined by 

This is appropriate because as mentioned above, {v ■ h)j = *v^. Applying 
Hodge star to both sides of equation (3.4) and replacing hj f everywhere 
we get the governing equations in terms of the flux and pressure, which can 
be written as 

k 

(3.8) / + -(*dop) = inM, 

A* 

(3.9) dn-if = (pto inM, 

(3.10) f = i'l ondM. 

Given k, fi, (j), ip and the boundary condition (3.10), the problem state- 
ment is to solve equations (3.8) and (3.9) for the flux / and pressure p. 
Equations (3.8)-(3.10) are the ones that we will discretize first in Section 4 
using the discrete operators defined in Section 2.4. An equivalent form for 
equation (3.8) obtained by applying Hodge star to both sides of (3.8) is 

(3.11) */ + (-l)"-^-dop = inM, 

and we will also discretize this equation to get an alternative formulation. 
Here the (—1)"^^ sign has come from the double application of Hodge star 
using equation (2.1). 
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4. Discretization of Equations 

Let K he a simplicial complex that approximates M and * K the circum- 
centric dual of K as defined in Section 2.2. Let L be the approximation of 
the boundary dM so that L consists of the (n — l)-dimcnsional boundary 
faces of K. The differential forms /, k, (f)u and V'7 in equations (3.8)- (3.10) 
are discretized as cochains and the operators d and * are replaced by their 
discrete counterparts described in Section 2.5. 

An important point to note when discretizing is the appropriate placement 
of the cochains. In particular we will place the discrete flux / on the (n — 
l)-dimensional primal simplices - edges in triangle mesh and triangles in 
tetrahedral mesh. Thus / G C"'~^(iC). From equation (3.8) this implies that 
the discrete version of * do p must also be placed on these primal simplices 
since kf/iis a scalar here. Thus * dop G C"^^(K) from which it follows that 
do p is a dual cochain and dop G C"~("~^)(*if), that is dop G C^{-kK) is a 
dual 1-cochain placed on the dual edges. This finally leads to the conclusion 
that discrete pressure is a dual 0-cochain, that is, p G C^{-kK) and thus 
the pressure must be placed at the circumcenters of the top dimensional 
simplices. 

Since /c is a dual 0-cochain wc must use the discrete operator dg to replace 
the exterior derivative in equation (3.8). Referring to the diagrams (2.6) 
and (2.7) it is clear that the discrete Hodge that should be used to replace 
* in equation (3.8) is (— l)"~^*^li, the (—1)""^ sign coming from expres- 
sion (2.5). Finally, since / G C"'~^{K) clearly the discrete exterior derivative 
d„_i will replace the smooth d„_i when discretizing equation (3.9). 

Thus the discretized equations corresponding to equations (3.8)- (3. 10) are 
the very similar looking 

(4.1) /+^((_i)-i*;^i^dSp)=0 ini^, 

(4.2) dn-if = (pu iuK, 

(4.3) / = V7 onL, 

where the unknowns and data are the cochains / G C^~^{K), p G C^{-kK), 
(f)u G C"^{K) and ■07 £ C"'~^{K) where the last one is carried by L. The 
matrix representation of equations (4.1) and (4.2) adjusted for the boundary 
condition (4.3) is the linear system to be solved which is described next. 
Let / be the vector representing the cochain / in the basis 

for C"-^^{K) described in Section 2.4. Recall that a* is the {n — l)-cochain 
that is on fij, the (n — 1) dimensional simplex number i. As mentioned 
in Section 2.4 the simplices are ordered, for example, in dictionary order 
[7]. Similarly, let k, (puj and ip'y be the vectors corresponding to the other 
quantities appearing in equations (4.1)-(4.3). 
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To obtain the linear system to solve we first write equations (4.1) and (4.2) 
in block matrix form using the matrix form of the operators and objects. 
This yields 



(4.4) 



(A;//x)(-l)("-i)M-_\[dS] 













' " 











where / is an AT^-i x iV„_i identity matrix, is an x zero matrix and 
[dg] is the matrix form of do- Then using the fact that [dg] = {—1)'^D^_-^ 
we get 



(4.5) 













' " 




p. 




0a; 



Assuming that the domain has only one connected component, the pres- 
sure is unique only up to a constant. Hence, pressure at a single point 
must be fixed to an arbitrary value to get a unique solution. To impose the 
boundary conditions the right hand side of equation (4.5) is adjusted for the 
known boundary fluxes. Such an adjustment is also done for the assumed 
pressure at one point. This is a standard procedure and it is described here 
briefly for completeness. The adjustments are done by taking the linear 
combination of the columns of the linear system matrix in (4.5) correspond- 
ing to the known / and k. The coefficients in the linear combination are the 
known / and k values. The result is subtracted from the right hand side. 
These columns and corresponding rows are then deleted from the matrices 
in equation (4.5). 

Equation (4.5) can be written in a simpler standard saddle point form 
[8] by starting from equation (3.11) instead of (3.8). This is equivalent to 
multiplying the first block row of (4.5) by — (/x/A;)M„_i and the resulting 
system is 

"-(Ai/A;)M„_i Dl_; 

(4.6) 



Dr. 









7" 




' " 








(pOJ 



To take into account the boundary conditions (4.3) and the assumed pressure 
at a point, the adjustments described above are applied to equation (4.6). 
The matrix on the left hand side of (4.6) is a saddle type matrix of the form 

'A B^^ 
B 

Even after the boundary fluxes and assumed pressure are taken into account 
the form of the matrix stays the same, using sub-matrices of A and B. 
Here the matrix ^ is a diagonal matrix with nonzero diagonal entries since 
A = {—fi/k)Mn-i and M„_i is the diagonal discrete Hodge matrix. One 
difference between our method and Raviart- Thomas finite elements is that 
in the latter case A is not a diagonal matrix. 
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The system obtained by adjusting equations (4.5) or (4.6) for boundary 
conditions and assumed pressure at a point can be solved for the remaining 
pressures and the unknown fluxes by using any of the standard techniques 
for saddle type matrices [8]. For the results in this paper we used the Schur 
complement reduction method [8] or a direct solve for the full matrix system. 
In our case Schur complement is particularly simple since the A matrix can 
be explicitly inverted trivially. 

4.1. Flux visualization. Once the flux has been determined, we visual- 
ize it by using Whitney interpolation as described in Section 2.6 to get a 
smooth (n — l)-form inside each n-simplex. We then obtain the correspond- 
ing velocity vector fields sampled at barycenters. The sampling could be 
done at any location or locations in the interior of the n-simpliccs, not just 
at the barycenter. Given a flux value, we can determine the velocity as fol- 
lows. The flux / is related to the velocity by / = *{v^) which implies that 
* / = = (-l)"-iv^ 

Consider first n = 2 and let the value of the Whitney interpolated flux 

1- form at a sampling point be adx + bdy where a and b are some constants. 
Then 

v" = — (* /) = — *{adx + b dy) = bdx — ady , 

which implies that for standard metric in E?, the velocity v at that point is 
the vector (6, —a). For n = 3 if the value of the Whitney interpolated flux 

2- form at a sampling point inside a tetrahedron is / = ady A dz + bdz A 
dx + cdx Ady then the associated velocity is given by 

v'" = * f = *{a dy A dz + b dz A dx + c dx A dy) = adx + bdy + cdz , 

which implies that for standard metric in M^, the velocity v is the vector 
(a, 6, c). 

5. Heterogeneous Permeability and Hodge Star 

In many physical problems the Hodge star (which is an operator depend- 
ing on the metric) appears as a material dependent operator. For example 
when Maxwell's equations are written in terms of differential forms, the elec- 
tric permittivity and the magnetic permeability are both Hodge stars [15]. 
The permittivity Hodge star relates the electric field 1-form to the electric 
induction 2-form and the permeability Hodge star relates the magnetic flux 
density 2-form to the magnetic field 1-form. 

In this section we rewrite the flux form of Darcy's law (equation (3.8)) in 
a form that permits its discretization when the permeability is a spatially 
dependent scalar quantity. The scalar permeability value is discretized as 
constant in each n-simplex. From these scalar values of the permeability, 
a spatially dependent discrete Hodge star operator is constructed. This in- 
volves combining the scalar permeabilities across an (n — l)-simplex in a 
weighted average that is suggested by DEC and described in this section. 
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The resulting matrix for this heterogeneous Hodge star is stiU a diagonal ma- 
trix. In this case, ensuring that the entries be positive requires a stronger 
condition than the Dclaunay condition for the homogenous Hodge star. For 
the heterogenous Hodge star we require that the n-simplices at the interface 
contain their circumcenters. Such a simplex is called a well-centered sim- 
plex. For the planar case well-centered triangles are acute angled. Methods 
for creating wcll-ccntcred meshes are discussed in [45, 46, 47]. For the het- 
erogenous Hodge star defined in this section we don't need the entire mesh 
to be well-centered, just the layer on either side of the interface, as long as 
the overall triangulation is Delaunay. This is much simpler than creation 
of well-centered meshes and not much harder than creation of Delaunay 
meshes. 

We start with the equations (3.8) - (3.10) with the exception that the 
term {k/ fj,){*dp) is replaced by (l//i)(**^ dp). Thus Darcy law (3.8) now 
becomes 

(5.1) f + -{*'' dp) =0 inM, 

where the Hodge star operator *^ is the spatially dependent Hodge star that 
depends on the permeability of the medium. The continuity equation and 
boundary condition stay the same as equations (3.9) and (3.10). This het- 
erogeneous Hodge star is discretized as a diagonal matrix where the diagonal 
entry corresponding to an (n — l)-simplex a"'~^ in the matrix {M^_i)~^ is 

|(t"-^| A;+|*cr"-in(T!f:| A;_|*(7"-^ncr!!:| 

(5.2) 



Here cr" and cr" are the two simplices that contain cr""^, and k+ and k- are 
the permeabilities in these n-simplices. The dual edge ★cj""^ points from 

o"_ into The expression |*(T"'~^ncr" | stands for the length of the portion 
of dual edge ★a""^ that lies in a" etc. For k^ = k- = k expression (5.2) 
reduces to 



k 



which is the corresponding diagonal entry of M^\, thus yielding the usual 
discretization of the homogeneous Hodge star scaled by k. 

The ratio on the right in expression (5.2) can be interpreted as a weighted 
average of permeabilities. Note that it is not a simple arithmetic or geomet- 
ric mean of the permeabilities. The weights are the same as the ones that 
were used for averaging piecewise constant vector fields along a shared face 
in [28]. In [28] it was shown that these are the unique weights that yield a 
discrete divergence theorem. See [28, Figure 5.4, Section 5.5 and Section 6.1] 
for more details. 
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6. Numerical Results 

We illustrate the performance of the proposed DEC based numerical 
method for Darcy flow using many standard test problems. In all the fig- 
ures here that show a velocity vector field, the flux / has been visualized 
as the corresponding vector field. As described in Section 4.1 the velocity 
vector field v is obtained from / by using the relationship = (—1)""^ * /. 
The velocity vector field is sampled at barycenters of top dimensional sim- 
plices and displayed as arrows based at barycenters. The pressure in most 
figures is displayed by plotting it against the x coordinate of the circumcen- 
ter which is where the pressures are defined. The only exceptions are the 
pressure plots in Figure 7 in which the pressure is displayed by coloring the 
triangle with a single color based on pressure value at the circumcenter. Fig- 
ures 4, 5, 7, 6, 11, and 12 were constructed from data generated by a Python 
implementation of our numerical method that used the PyDEC software [7] ; 
and Figures 8, 9 and 10, were generated by a MATLAB implementation of 
our numerical method. 

6.1. Patch tests in 2D. The first results shown in Figure 4 are for patch 
tests [32]. It is desirable that for simple meshes a numerical method for 
Darcy flow should reproduce constant velocity and linear pressure exactly 
up to machine precision. The boundary condition in these tests is derived 
from a constant horizontal velocity (1,0). Thus, for example, in the square 
domain v ■ n = ip = —1 oia. the left edge, ip = 1 on the right edge, and on 
the top and bottom edges of the square. Keeping in mind the orientations 
and the definition of 7, when the discretized equations (4.1)-(4.3) are used, 
/ = 1 on the left and right edges and on the top and bottom edges of the 
square. 

We now explain the sign of / in more detail. Suppose the bottom left and 
top left corner vertices of the square domain in Figure 4 are labeled vq and 
v\ and the edge between them is oriented from vq to v\. We will use the 
name a for this oriented edge [wq, v\] and denote the vector from vq to v\ by 
a. Assume also that the square, and hence the triangle to which a belongs, 
is oriented counterclockwise. We use the same name / for the 1-cochain 
/ of equations (4.1)-(4.3) and the 1-form / of equations (3.8)-(3.10), but 
to be more precise the cochain / should be referred to as R{f) where R is 
the de Rham map of Section 2.4. The following calculation explains why 
{R{f),a) = +1 in this setting. 

{R{f),a) = / Vt = {v ■ h)-f = {v ■ n) / 7 

= (-1) / 7 = -7(<?) = a) = +1 . 

Here the third equality follows from the fact that v ■ h is constant along 
a, the fifth equality is true because a is a straight line, and uj{h,a) = — 1 
because the length of a is 1 and the basis (n, a) for the plane is oriented 
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clockwise, which is opposite of the orientation of the square. If a had been 
oriented from v\ to vo instead, the value of {R{f),a) would have been — 1 
instead. 

In this test, we are also given that = 0, so that divv = (equivalently, 
d / = 0) in the domain. The parameters k and /j, are 1, and g = so there is 
no external forcing. This example is constructed by starting with pressure 
p = —x + c for some constant c. Then it follows that v = (1,0) everywhere 
and div f = as given. The numerical method is given the (p and ^ and the 
pressure and flux is computed using the method. Although such patch tests 
do not guarantee that a method is high quality [6] , it is a convenient way to 
find problems with methods, as shown in Figure 1. If a method fails such a 
simple test, it is probably unsuitable for the problem. 

The relative errors for the nonzero pressures arc less than 3 x 10^^® for 
the square and less than 7 x 10^^^ for the hexagon shown in Figure 4. Thus 
for these simple meshes the relative error is close to machine precision. The 
same test is repeated for a larger mesh of a square domain in Figure 5 for 
which the relative error in pressure is less than 9 X 10" Fi gure 6 shows the 
result of same patch test on a random Delaunay mesh, in which 37 vertices 
were placed randomly, with no regard to mesh quality. 

6.2. Known solution vi^ith nonzero source term. In the next test we 
compare the solution computed using our method with an analytically known 
solution. Again the parameters k and are 1 and g = 0. The solution is 
constructed by starting with pressure p = cos(7rx) cos(7ry) from which an 
expression for v is derived. From v one computes the divergence to derive 
that the source/sink term is ^ = 27r^ cos(7ra;) cos(7r?/). The boundary data 

is constructed from the analytically computed v. The numerical method 
is given (p and tp and used to compute v and k in the domain from that. 
Figure 7 shows a comparison of the computed pressure and velocity with 
the analytical solution. 

6.3. Discontinuous permeability. One of the new results of our approach 
is a Hodge star that allows the discretization of the equations in the case 
of spatially varying scalar permeability. The permeability is taken to be 
constant in each n-simplex. The resulting heterogeneous diagonal Hodge 
star matrix was defined in Section 5. An important aspect of any numerical 
method for Darcy flow is how well it can address discontinuities in perme- 
ability. Our heterogeneous discrete Hodge star allows us to test this aspect 
of our method. Figure 8 shows the results of a patch test with discontinuous 
permeability in two adjoining domains. The domain is a rectangle in which 
the triangles in the left half are given a permeability of ki and those in the 
right half are given a permeability of k2- The fluid flows in from the left 
and exits from the right. The pressure should be a piecewise linear function 
whose slope depends on the permeability and the velocity in the domain 
should be constant. This is demonstrated in the results from our method. 
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6.4. Layered medium. Another common test in the Darcy flow literature 
is when the discontinuities in the permeabihty vary across the flow rather 
than along it. Such a medium is typically called a layered medium in which 
the various layers are given different permeabilities. In our layered medium 
computation we perform two tests. In one the permeabilities alternate be- 
tween 1 and 5 and in another they alternate between 1 and 10. The fluid 
comes into the domain from the left and exits from the right as in our patch 
tests. The velocity should be horizontal and constant in a layer. It should 
be larger in the low permeability layers as computed by our method and 
shown in Figure 9. The correct pressure profile should be a linear function 
of X and this is seen in Figure 10. 

6.5. Patch tests in 3D. Many numerical formulations perform well in 2D, 
but their natural extensions to 3D do not perform well. Herein we show 
that the proposed formulation performs well even in 3D. To illustrate this 
we considered two different computational domains, which are shown in 
Figures 11 and 12. The first domain is a polyhedron with 16 tetrahedra and 
the second domain is a cube with 244 tetrahedra. The analytical solution 
is constant velocity along x direction, and pressure linearly varying along x 
direction. The obtained numerical results are plotted in Figures 11 and 12 
which shows that the numerical method performed well. For example, the 
latter figure shows that the relative error in pressure is less than 2 x 10~^^. 

6.6. Circumcenter versus barycenter. One of the main focuses of this 
paper is on structure preserving numerical methods, generated by a system- 
atic application of DEC. The proposed method has both local and global 
mass balance properties. In addition, we show that the proposed method 
can also exactly represent linear variation of pressure within a domain. We 
also highlight the need for careful choice of locations for the pressure to get 
better numerical solutions. DEC naturally provides the locations of pres- 
sure that conserve local and global mass balance, and also exactly represent 
linear variation of pressure. In DEC, for each element, pressure is located 
at circumcenter. The fact that each (n — l)-simplex and its circumcentric 
dual edge are orthogonal is important here. One common choice used in 
the literature as location points for pressure are barycenters. In Figure 3 
we show that choice of barycenter cannot exactly represent linear variation 
of pressure (along with local and global mass balance) whereas location of 
pressure at circumcenter can. This also illustrates how a geometrical view of 
PDEs can provide good insights into developing novel structure preserving 
numerical methods. 

7. Conclusions and Future Work 

Our goal in this paper has been to introduce a numerical method for Darcy 
flow derived using discrete exterior calculus (DEC). We have shown that this 



20 



A. N. HIRANI, K. B. NAKSHATRALA, AND J. H. CHAUDHRY 



0.5 




(D 








0.2 



0.4 



0.6 



0.8 



X 



Figure 3. Locating pressures at barycenters does not work. 
Shown here is the pressure computed for a square domain in which 
the boundary condition is given by constant horizontal velocity 
pointing to right. Thus velocity in the domain should be constant 
and pressure linear which is reproduced when pressures are located 
at circumcenters as dictated by DEC. 



approach results in a unified derivation of methods for both two- and three- 
dimensions. We have also demonstrated the numerical performance of these 
methods. For example, the method is shown to pass several patch and other 
standard test problems in both two and three dimensions. Our DEC based 
approach is intrinsic, in the sense that it involves quantities and operations 
that do not depend on how the mesh is embedded in M^. As a result, it 
should be easy to use the method for Darcy flow on a curved surface using 
a discretization of such a surface. This is one avenue for future research. 

Even in the case of domains that are open subsets of or the re- 
quired computations such as circumcenter calculation, and computation of 
volumes, areas and lengths, can be built from linear algebra operations which 
simplifies the implementation as demonstrated in [7]. 

The philosophy of DEC is to discretize the operators and objects in such 
a way that the their properties in smooth calculus have discrete analogs. As 
a result, for example, in our method mass balance is satisfied both locally 
(that is, in each element) and globally because of the way the fluxes are 
represented and because Stokes' theorem is true by definition in DEC. 

A pedagogical advantage of the DEC approach is that once the language 
of exterior calculus has been mastered and the discretizations understood, 
the translation from a PDE to its corresponding discrete equation, and from 
there to the matrix form are trivial steps. Compare for example the smooth 
equation (3.8), the corresponding discrete equation (4.1) and the first block 
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row of the matrix equation (4.5). There is also a clean separation in the ma- 
trix form equation (4.6) between the parts depending only on the mesh con- 
nectivity (the top right and bottom left blocks) and that depending on how 
the physical variable are related (the top left part) which usually involves 
constitutive relationships as well. These depend on material measurements, 
and hence the inaccuracies of measurement of material properties does not 
corrupt that part of the matrix that is topological. In this aspect our method 
shares this good property of some mixed finite element method formulations. 

In addition to applying DEC to Darcy flow, in this paper we have also 
developed a discretization of Hodge star for non-homogcncous medium. Wc 
used this in examples in which the permeability varies across the mesh, either 
along the flow or transversal to the flow. This discretization of a spatially 
varying Hodge star should be useful in other PDEs as well. 

There are many avenues for future research even within Darcy flow. Flow 
on curved surfaces has already been mentioned before. Another direction 
for further work is the proper DEC discretization of anisotropic permeabil- 
ity. Finally the analytical treatment of convergence and stability properties 
remains to be done. For this, it might be useful to formulate this method 
as a finite element method. 
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Figure 4. Results from patch tests. The boundary conditions are 
obtained from a velocity of 1 in positive x direction. Thus, for the 
square domain, in equation (3.3), ^ = — 1 on left edge, 1 on right 
edge, and on top and bottom edges. For both the square and 
hexagon domains it is also given is that (/) = in (3.2), constants 
A: = 1, /i = 1 and external body force acceleration g = in (3.1). 
The pressure should be linear as shown in the right figures. The 
left figures show that the velocity is constant, as expected. The 
velocity displayed is interpolated from the flux through each edge, 
using Whitney 1-form and sampled and converted to a vector. 
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Figure 5. Patch test with larger mesh and with the same bound- 
ary conditions, parameters and other data as the square in Fig- 
ure 4. 
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Figure 6. Patch test using a Delaunay triangulation of randomly 
placed 37 vertices in a square. Several triangles are obtuse angled 
and no attempt has been made to make the mesh have high quality. 
The boundary conditions, parameters and other data are same as 
for the square in Figure 4. 
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Figure 7. Comparison of analytically known solution with a solu- 
tion computed using the method developed here. The top left fig- 
ure shows the pressure part of the analytical solution to the Darcy 
law problem with k — 1, fi — 1, g — 0, and (f) — liP' cos(7ra;) cos(7r?/). 
This implies p = co%(txx) cos^ny) which is shown overlaid on the 
mesh. In the left bottom figure we show the pressure computed 
with the numerical method proposed here. Each triangle is colored 
by the pressure at its circumcenter. The top right figure shows the 
analytically computed vector field, and the numerically computed 
vector field is shown in the right bottom figure, which is obtained 
from the flux by interpolation of 1-forms. 
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Figure 8. In this test the left and right halves of the mesh are 
given different permeabilities, fci on the left and fc2 on right. Thus 
the permeability jumps across the middle vertical edge at x = 0.5. 
The values that we use for (fci,A:2) are (1,1), (1,2), (1,10) and 
(1, 100). The boundary condition is derived from velocity (1, 0) as 
in Figure 4. The top figure shows the computed velocity interpo- 
lated from flux for fci = 1 and k2 = 10. All other value pairs also 
result in a constant velocity solution as expected. The pressure 
is piecewise linear with a jump at the discontinuity at a; = 0.5 as 
shown in the bottom figure. 
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Figure 9. Layered medium with 2 different permeability patterns. 
The domain has 5 layers with alternating permeability. In the top 
figure the permeability fc, from bottom layer to top is 5, 10, 5, 10 
and 5. In the bottom figure the permeability k is 1, 10, 1, 10 and 
1. The computed flux is visualized as a vector field. 
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Figure 10. Pressure for the layered medium show in Figure 9. 
The pressure is hnear as expected. The top and bottom figures 
correspond to the top and bottom figures in Figure 9. 
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Figure 11. Patch test in 3D. The mesh shown has 16 tetrahedra. 
The fluid velocity on the boundary is from negative x to positive 
X direction, with no y or z components. The fluxes (2-cochains) 
across the internal faces are computed using our method and in- 
terpolated using the Whitney map. This is then sampled at the 
circumcenters, converted into a vector field and plotted as arrows. 
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Figure 12. Patch test on a cube. The mesh used for this cube has 
244 tetrahedra. For clarity, the tetrahedra in the cube have not 
been shown. The fluid velocity on the boundary is from negative x 
to positive x direction, with no y or z components. The fluxes (2- 
cochains) across the internal faces are computed using our method 
and interpolated using the Whitney map. This is then sampled 
at the circumcenters, converted into a vector field and plotted as 
arrows. In the relative error plot for pressure, the data with 
exact pressure has been removed. 



